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ABSTRACT 

We simulate the growth of galaxies and their central supermassive black holes by 
implementing a suite of semi-analytic models on the output of the Millennium Run, 
a very large simulation of the concordance ACDM cosmogony. Our procedures fol- 
low the detailed assembly history of each object and are able to track the evolution 
of all galaxies more massive than the Small Magellanic Cloud throughout a volume 
comparable to that of large modern redshift surveys. In this first paper we supple- 
ment previous treatments of the growth and activity of central black holes with a new 
model for 'radio' feedback from those AGN that lie at the centre of a quasistatic X-ray 
emitting atmosphere in a galaxy group or cluster. We show that for energetically and 
observationally plausible parameters such a model can simultaneously explain: (i) the 
low observed mass drop-out rate in cooling flows; (ii) the exponential cut-off at the 
bright end of the galaxy luminosity function; and (iii) the fact that the most massive 
galaxies tend to be bulge-dominated systems in clusters and to contain systematically 
older stars than lower mass galaxies. This success occurs because static hot atmo- 
spheres form only in the most massive structures, and radio feedback (in contrast, for 
example, to supernova or starburst feedback) can suppress further cooling and star 
formation without itself requiring star formation. We discuss possible physical mod- 
els which might explain the accretion rate scalings required for our phenomenological 
'radio mode' model to be successful. 

Key words: cosmology: theory, galaxies: formation, galaxies: evolution, galaxies: 
active, galaxies: cooling flows, black hole physics 



1 INTRODUCTION 

The remarkable agreement between recent measurements of 
cosmic structure over a wide range of length- and time- 
scales has established a standard paradigm for structure 
formation, the ACDM cosmogony. This model can simulta- 
neously match the microwave background fluctuations seen 
at z ~ 1000 (e.g. ^Dcrgcl ot al. 2003), the p ower spectrum 
of th e low redshift galaxy distribution (e.g. iPercival et alJ 
l2002l: iTeemark et"al] |2004) . the nonlinear mass distribu- 
tion at low redshift as characterised by cosmic shear (e.g. 
Van Waerbeke et a l. 2002j and the structure seen in the 
z = 3 Ly Q forest (e.g. Mandelbau m et al. 2003). It also re- 
produces the present acceleration of t he cosmic expansion in- 
ferred from supe rnova observations JPerlmutter et alJll999l: 
iRiess et al.lll99^ . and it is consistent with the mass budget 
inferred for the present universe from the dynamics of large- 



scale structure ("Peacock et al.'"200l|), the baryon fraction in 
rich clusters (|Whit^^t^j^93) and the theory of Big Bang 
nucleosynthesis""]o]^^^^]|2000). A working model for the 
growth of all structure thus appears well established. 

In this cosmogony, galaxies form when gas condenses at 
the centres of a hierarchically merging populat i on of dark 
haloes, as originally proposed bv I White fc Ree^ lll97d) . At- 
tempts to understand this process in detail have consistently 
run into problems stemming from a mismatch in shape be- 
tween the predicted distribution of dark halo masses and 
the observed distribution of galaxy luminosities. Most stars 
are in galaxies of Milky Way brightness; the galaxy abun- 
dance declines exponentially at brighter luminosities and in- 
creases sufficiently slowly at fainter luminosities that rela- 
tively few stars are in dwarfs. In contrast, the theory pre- 
dicts a much broader halo mass distribution - a constant 
mass-to-light ratio would produce more high and low lu- 
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minosity galaxies than are observed while underpredicting 
the number of galaxies like the Milky Way. Attempts to 
solve these problems initially invoked cooling inefficiencies 
to reduce gas condensation in massive systems, and super- 
nova fee dback to reduce star formation efficiency in low mass 
systems JWhite fc ReesllOTStlWhite fc Frenkll99lfl . Forma- 
tion of dwarfs may also be suppressed by photoionisa tion 
heating jEfstathiou 'l992Y As'Th oul fc Weinberd ^ll995^ em- 
phasised, cooling effects alone are too weak to produce the 
bright end cut-off of the luminosity function, and recent at- 
tempts to fit observed luminosity functions have been forced 
to in clude additional fee dback processes in massive systems 
(e.g. iBenson et aP '2003V In this paper we argue that ra- 
dio sources may provide the required feedback while at the 
same time providing a solution to two other long-standing 
puzzles. 

An important unanswered question is why the gas at 
the centre of most galaxy clusters is apparently not con- 
densing and turning into stars when the observed X-ray 
emission implies a cooling time which is much shorter than 
the age of the system. This cooling flow puzzle was noted 
as so on as the first X-ray maps of clusters becam e avail- 
able JCowie fc Binnevlll977l: iFabian fc Nulsenlll977l) and it 
was made more acute when X-ray spectroscopy demon- 
strated that very little gas is cooling through temperatures 
even a factor of three below that of th e bulk of the gas 
JPeterson et alJl200ll: iTamura et"ai]l200l[). A clue to the so- 
lution may come from the observation IIBurns et al1ll98lll 
that every cluster with a strong cooling flow also contains 
a mas sive and active central radio galaxy. Tabor & Binncv 
il993 h suggested that radio galaxies might regulate cool- 
ing flows, and this idea has gained considerable recent 
support from X-ray maps which show direct evidence for 
an interaction between radio lobes and the intracluster 
gas (_Fabian ct al. 2003; McNamara et al. 2Q00, 2005). A 
number of authors have suggested ways in which the ra- 



dio source might replace the thermal energy lost 


to X- 


ray emission ijBinnev & Tabodll995l IChurazov et al.l 


2002 


Briieeen & Kaiser' '200^; 'Ruszkowski fc BeeelmanI 


2002 


Kaiser & Binncv 2003; Omnia ct al. 2004)- We do not focus 



on this aspect of the problem here, but rather demonstrate 
that if the scaling properties of radio source feedback are set 
so they can plausibly solve the cooling flow problem they 
induce a cut-off at the bright end of the galaxy luminosity 
function which agrees well with observation. 

Another puzzling aspect of the galaxy population is 
the fact that the most massive galaxies, typically ellipti- 
cals in clusters, are made of the oldest stars and so fin- 
ishe d their star formation earlier than lower mass galax- 
ies (iBender fc Saglial ll999D . Confirming evidence for this 
comes from look-back studies which show that both star- 
formation and AGN activity take place more vigorously 
and in higher mass objec ts at redshifts of 1 to 2 than in 
the p r esent Univer s e (e.g . IShaver et al.lll99^ iMadau et alJ 
119961) . ICowie et alJ il996l) termed this phenomenon 'down- 
sizing', and prima facie it conflicts with hierarchical growth 
of structure in a ACDM cosmogony where massive dark 
halo es assemble at lowe r redshift than lower mass haloes 
(e.g. iLacev fc Coi3ll993l) . This puzzle is related to the pre- 
vious two; the late-forming high mass haloes in ACDM corre- 
spond to groups and clusters of galaxies, and simple theories 
predict that at late times their central galaxies should grow 



to masses larger than those observed through accretion from 
cooling flows. In the model we present below, radio galaxies 
prevent significant accretion, thus limiting the mass of the 
central galaxies and preventing them from forming stars at 
late times when their mass and morphology can still change 
through mergers. The result is a galaxy luminosity function 
with a sharper high-mass cut-off in which the most massive 
systems are red, dead and elliptical. 

To make quantitative predictions for the galaxy popu- 
lation in a ACDM universe it is necessary to carry out sim- 
ulations. Present numerical capabilities allow reliable sim- 
ulation of the coupled nonlinear evolution of dark mat- 
ter and diffuse gas, at least on the scales which determine 
the global properties of galaxies. Once gas cools and con- 
denses into halo cores, however, both its structure and the 
rates at which it turns into stars and feeds black holes are 
determined by small-scale 'interstellar medium' processes 
which are not resolved. These are usually treated through 
semi-analytic recipes, parameterised formulae which encap- 
sulate 'subgrid' physics in terms of star formation thresh- 
olds, Schmidt 'laws' for star formation, Bondi models for 
black hole feeding, etc. The form and the parameters of 
these recipes are chosen to reproduce the observed sys- 
tema tics of star form ation and AGN activity in galaxies 
(e.g. lKennicut3 Il998() . With a well-constructed scheme it 
is possible to produce stable and numerically converged 
simulation s which mimic real star-formi ng galaxies remark- 
ably well JSprineel fc Hernauisdl2003a^ . In strongly star- 
forming galaxies, massive stars and supernovae produce 
winds which redistri bute energy, ma s s and heavy element s 
over large regions iHeckman et alJ Il990l : iMartinI 1199911 . 
Even stronger feedbac k is possible, in principle, from AGN 
jBegelman et al.lll99iri . In both cases the determining pro- 
cesses occur on very small scales and so have to be included 
in simulations through parametrised semi-analytic models. 
Unfortunately, the properties of simulated galaxies turn out 
to depend strongly on how these unresolved star-formation 
and feedback processes are treated. 

Since the diffuse gas distribution and its cooling onto 
galaxies are strongly affected by the description adopted for 
the subgrid physics, every modification of a semi-analytic 
model (or of its parameters) requires a simulation to be re- 
peated. This makes parameter studies or tests of, say, the ef- 
fects of different AGN feedback models into a very expensive 
computing exercise. A cost-effective alternative is to repre- 
sent the behaviour of the diffuse gas also by a semi-analytic 
recipe. Since the dark matter couples to the baryons only 
through gravity, its distribution on scales of galaxy haloes 
and above is only weakly affected by the details of galaxy for- 
mation. Its evolution can therefore be simulated once, and 
the evolution of the baryonic component can be included 
in post-processing by applying semi-analytic models to the 
stored histories of all dark matter objects. Since the sec- 
ond step is computationally cheap, available resources can 
be used to carry out the best possible dark matter simula- 
tion, and then many parameter studies or tests of alterna- 
tive models can be carried out in post-processing. T his sim- 
ulatio n approach was first introduced by tKauffmann et alJ 
il99il) and it is the approach we apply here to the Mil- 
lennium Run, the largest calculation to date of the evo- 
lution of structure in the concordance ACDM cosmogony 
JSprineel et alJl2005b^ . 
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This paper is organised as follows. In Section |5| we de- 
scribe the Millennium Run and the post-processing we car- 
ried out to construct merging history trees for all the dark 
haloes within it. Sectionj^presents the model for the forma- 
tion and evolution of galaxies and their central supermas- 
sive black holes that we implement on these merging trees. 
Section 2] describes the main results of our modelling, con- 
centrating on the influence of 'radio mode' feedback on the 
properties of the massive galaxy population. In Section |S] 
we discuss physical models for black hole accretion which 
may explain the phenomenology required for our model to 
be successful. Finally, Section |n] summarises our conclusions 
and suggests some possible directions for future investiga- 
tion. 



2 THE DARK MATTER SKELETON: THE 
MILLENNIUM RUN 

Our model for the formation and evolution of galaxies and 
their central supermassive black holes is implemented on 
top of the Millennium Run, a very large dark matter simu- 
lation of the concordance ACDM cosmology with 2160^ ~ 
1.0078 X 10^° particles in a periodi c box of 500 h~ ^ Mpc o n 
a side. A full description is given in ISprineel et ai]ll2005bD : 
here we summarise the main simulation characteristics and 
the definition and construction of the dark matter merging 
history trees we use in our galaxy formation modelling. The 
dark matter distribution is illustrated in the top panel of 
Fig. □ for a 330 x 280 x 15/i"^Mpc slice cut from the full 
volume. The projection is colour coded by density and lo- 
cal velocity dispersion, and illustrates the richness of dark 
matter structure for comparison with structure in the light 
distribution to which we will come later. Dar k matter plots 
on a w ider range of scales may be found in ISpringel et alJ 
J2005d) . 

2.1 Simulation characteristics 

We adopt cosmological parameter val ues consistent with 
a combined analysis of th e 2dFGRS JColless et alJl200lll 
and first-year WMAP data JSoergel et al.l2003l : ISeliak et alJ 
^005). They are = i^dm + = 0.25, = 0.045, 
h = 0.73, = 0.75, n — 1, and erg = 0.9. Here Q^n de- 
notes the total matter density in units of the critical density 
for closure, pait = 3iii'o/(87rG). Similarly, f2b and Ha denote 
the densities of baryons and dark energy at the present day. 
The Hubble constant is given as Ho = 100/ikms"^Mpc"\ 
while erg is the rms linear mass fluctuation within a sphere 
of radius 8/i~^Mpc extrapolated to 2 = 0. 

The chosen simulation volume is a periodic box of 
size 500/i~^Mpc, which implies a particle mass of 8.6 x 
10*/i~^Mq. This volume is large enough to include inter- 
esting objects of low space density, such as quasars or rich 
galaxy clusters, the largest of which contain about 3 mil- 
lion simulation particles at z — 0. At the same time, the 
mass resolution is sufficient that haloes that host galaxies 
as faint as 0.1 are typically resolved with at least ~ 100 
particles. Note that although discreteness noise significantly 
affects the merger histories of such low mass objects, the 
galaxies that reside in halos with ^ 100 particles are usu- 
ally sufficiently far down the luminosity function that any 



uncertainty in their properties has little impact on our re- 
sults or conclusions. 

The initial conditions at z = 127 were created by dis- 
placing particles from a homogeneous, 'glass-like' distribu- 
tion (White 1 99^) using a Gaussian random field with a 
ACDM linear po wer spectrum as given by t he Boltzmann 
code CMBFAST JSeliak fc Zaldarriagal [19961) . The simula- 
tion was then evolved to the present epoch using a leapfrog 
integration scheme with individual and adaptive time steps, 
with up to 11 000 time steps for individual particles. 

The simulation itself was carr ied out with a spe- 
cial version of the GADGET-2 code JSpringel et alj|2001bt 
ISpringeJiiool optimised for very low memory consump- 
tion so that it would fit into the nearly 1 TB of physi- 
cally distributed memory available on the parallel IBM p690 
computer^ used for the calculation. The computational al- 
gorithm used the 'TreePM' method (.Xu 1995; Bode ct aT] 
l2000tlBaglall2002^ to evaluate gravitational forces, combin- 
ing a hierarchical mu ltipole expansion, or 'tree' algorithm 
iBarnes fc Hud Il986l) , and a classical, Fourier t ransform 
particle-mesh method jHocknev fc EastwoodlllOS j) . An ex- 
plicit force-split in Fourier-space produces a very nearly 
isotropic force law with negligible force errors at the force 
matching scale. The short-range gravitational force law is 
softened on comoving scale 5 /i~^kpc (Plummer-equivalent) 
which may be taken as the spatial resolution limit of the 
calculation, thus achieving a dynamic range of 10^ in 3D. 
The calculation, performed in parallel on 512 processors, 
required slightly less than 350 000 processor hours of CPU 
time, or 28 days of wall-clock time. 



2.2 Haloes, substructure, and merger tree 
construction 

Our primary application of the Millennium Run in this pa- 
per uses finely resolved hierarchical merging trees which en- 
code the full formation history of tens of millions of haloes 
and the subhaloes that survive within them. These merging 
history trees are the backbone of the model that we imple- 
ment in post-processing in order to simulate the wide range 
of baryonic processes that are important during the forma- 
tion and evolution of galaxies and their central supermassive 
black holes. 

We store the full particle data between z = 20 and z = 
at 60 output times spaced in expansion factor according to 
the formula 

log (1 + 2„) = n(n -f 35)/4200 . (1) 

This spacing is 'locally' logarithmic but becomes smoothly 
finer at lower redshift, with a temporal resolution by red- 
shift zero of approximately 300 Myr. Additional outputs are 
added at z — 30, 50, 80, 127 to produce a total of 64 snap- 
shots in all. We note that each snapshot has a total size in 
excess of 300 GB, giving a raw data volume of nearly 20 TB. 

Together with each particle coordinate dump, the sim- 
ulation code directly produces a friends-of-friends (FOF) 
group catalogue on the fiy and in parallel. FOF groups are 



^ This computer is operated by the Computing Centre of the 
Max-Planck Society in Garching, Germany. 
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125 MpC/h 





the Millennium Run. For the dark matter distribution, intensity encodes surface density and colour encodes local velocity dispersion. 
For the light distribution, intensity encodes surface brightness and colour encodes mean B— V colour. The linear scale is shown by the 
bar in the top panel. 
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defined as equivalence classes where any pair of two parti- 
cles is placed into the same group if their mutu al separation 
is less than 0.2 of the mean particle separation to a vis et alJ 
0,98^. This criterion combines particles into groups with 
a mean overdensity of about 200, corresponding approxi- 
mately to that expected for a virialised group. The group 
catalogue saved to disk for each output only kept groups 
with at least 20 particles. 

High-resolution simulations like the present one ex- 
hibit a rich substructure of gravitationally bound dark mat- 
ter subhaloes orbit ing within larger virialised halos (e.g. 
iGhigna et allllQQ^ . The FOF group-finder built into our 
simulation code identifies the haloes but not their subha- 
los. Since we wish to follow the fate of infalling galaxies and 
their halos, and these are typically identifiable for a sub- 
stantial time as a dark matter subhalo within a FOF halo, 
we apply in post-processing an imp roved and extended ver - 
sion of the SUBFIND algorithm of ISoringel et al.1 ll2001all . 
This computes an adaptively smoothed dark matter density 
field within each halo using a kernel-interpolation technique, 
and then exploits the topological connectivity of excursion 
sets above a density threshold to identify substructure can- 
didates. Each substructure candidate is subjected to a grav- 
itational unbinding procedure. If the remaining bound part 
has more than 20 particles, the subhalo is kept for further 
analysis and some of its basic physical properties are deter- 
mined (angular momentum, maximum of its rotation curve, 
velocity dispersion, etc.). After all subhaloes are identified 
they are extracted from the FOF halo so that the remain- 
ing featureless 'background' halo can also be subjected to 
the unbinding procedure. This technique, however, neglects 
the fact that substructures embedded within a halo help to 
bind its material, and thus removing them can, in princi- 
pal, unbind some of the FOF halo that may otherwise still 
be loosely bound. We accept this small effect for techni- 
cal reasons related to the robustness of our halo definition 
procedures and the need for unambigous particle-subhalo 
assignments in our data structures. We note that the to- 
tal mass in substructures is typically below 10% and often 
substantially smaller. Thus any bias in the bound mass of 
the parent halo due to additional unbinding is always very 
small. 

To compute a virial mass estimate for each FOF halo we 
use the spherical-overdensity approach, where the centre is 
determined using the minimum of the gravitational potential 
within the group and we define the boundary at the radius 
which encloses a mean overdensity of 200 times the critical 
value. The virial mass, radius and circular velocity of a halo 
at redshift z are then related by 



Mvir 



10GH{z) 



(2) 



where H{z) is the Hubble constant at redshift z. 

At 2 = our procedures identify 17.7 x 10'' FOF groups, 
down from a maximum of 19.8 X lO" at 2 = 1.4 when groups 
were more abundant but of lower mass on average. The z — Q 
groups contain a total of 18.2 x 10® subhaloes, with the 
largest FOF group containing 2328 of them. (Note that with 
our definitions, all FOF groups contain at least one subhalo, 
the main subhalo which is left over after removal of any 
substructure and the unbound component. We require this 
main subhalo to contain at least 20 particles.) 



Having found all haloes and subhaloes at all output 
snapshots, we then characterise structural evolution by 
building merging trees that describe in detail how haloes 
grow as the universe evolves. Because structures merge hier- 
archically in CDM universes there can be several progenitors 
for any given halo, but in general there is only one descen- 
dant. (Typically the cores of virialised dark matter struc- 
tures do not split into two or more objects.) We therefore 
construct merger trees based on defining a unique descen- 
dant for each halo and subhalo. This is, in fact, sufficient to 
define the entire merger tree, since the progenitor informa- 
tion then follows impl icitly. Further details can be found in 
ISprineel et"ai] (l2005b^ . 

We store the resulting merging histories tree by tree. 
Since each tree contains the full formation history of some 
z — Q halo, the physical model for galaxy formation can be 
computed sequentially tree by tree. It is thus unnecessary 
to load all the history information on the simulation into 
computer memory at the same time. Actually, this would 
be currently impossible, since the combined trees contain a 
total of around 800 million subhaloes for each of which a 
number of attributes are stored. Note that, although evolv- 
ing the galaxy population sequentially in this way allows us 
to consider, in principle, all interactions between galaxies 
that end up in the same present-day FOF halo, it does not 
allow us to model longer range interactions which might take 
place between galaxies that end up in different FOF halos 
at z — Q. 



3 BUILDING GALAXIES: THE 
SEMI-ANALYTIC MODEL 

3.1 Overview 

In the following sub-sections we describe the baryonic 
physics of one particular model for the formation and evo- 
lution of galaxies and of their central supermassive black 
holes. A major advantage of our simulation strategy is that 
the effects of parameter variations within this model (or 
indeed alternative assumptions for some of the processes) 
can be explored at relatively little computational expense 
since the model operates on the stored database of merger 
trees; the simulation itself and the earlier stages of post- 
processing do not need to be repeated. We have, in fact, 
explored a wide model and parameter space to identify our 
current best model. We summarise the main parameters of 
this model, their 'best' values, and their plausible ranges in 
Table Q These choices produce a galaxy population which 
matches quite closely many observed quantities. In this pa- 
per we discuss the field galaxy luminosity-colour distribu- 
tion; the mean stellar mass - stellar age relation; the TuUy- 
Fisher relation, cold gas fractions and gas-phase metallicities 
of Sb/c spirals; the colour - magnitude relation of ellipticals; 
the bulge mass - black hole mass relation; and the volume- 
averaged cosmi£sta£jornm^^ and black hole accretion his- 
tories. In lSpringel et alj i2005lJ) we also presented results for 
galaxy correlations as a function of absolute magnitude and 
colour, for the baryonic 'wiggles' in the large-scale power 
spectrum of galaxies, and for the abundance, origin and 
fate of high redshift supermassive black holes which might 
correspond to t he 2 ~ 6 quasars discovered by the SDSS 
jFan et al.ll200lh 
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Table 1. A summary of our main model parameters and their best values and plausible ranges, as described in the text. Once set, these 
values are kept fixed for all results presented in this paper, in particular for models in which AGN feedback is switched off. 



parameter 


description 


best value 


plausible range 


ft 

Zq, Zr 


cosmic baryon fraction 
redshift of reionization (|^^^ 


0.17 

O, ( 


fixed 
fixed 


/bh 
«;agn 


merger cold gas BH accretion fraction fi|c{.4.1i 
quiescent hot gas BH accretion rate (MQyr^^) (4'iA.'2t 


0.03 

6 X lO-'^ 


0.02 - 0.04 
(4 - 8) X 10-6 


«SF 


star formation efficiency (4'i.5t 


0.07 


0.05 - 0.15 


Cdisk 
^halo 

7oj 


SN feedback disk reheating efficiency ('83. 6i 
SN feedback halo ejection efficiency 
ejected gas reincorporation efficiency ('i|3.(il 


3.5 
0.35 
0.5 


1-5 
0.1 - 0.5 
0.1 - 1.0 




major merger mass ratio threshold (t'i.Jt 


0.3 


0.2 - 0.4 


R 
Y 


instantaneous recycled fraction of SF to the cold disk f il3.9i 
yield of metals produced per unit SF f i|3.9l 


0.3 
0.03 


0.2 - 0.4 
0.02 - 0.04 



In our model we aim to motivate each aspect of the 
physics of galaxy formation using the best available obser- 
vations and simulations. Our parameters have been chosen 
to reproduce local galaxy properties and are stable in the 
sense that none of our results is critically dependent on any 
single parameter choice; plausible changes in one parame- 
ter or recipe can usually be accommodated through adjust- 
ment of the remaining parameters within their own plausible 
range. The particular model we present is thus not unique. 
Importantly, our model for radio galaxy heating in cooling 
flows, which is the main focus of this paper, is only weakly 
influenced by the remaining galaxy formation and black hole 
growth physics. This is because our radio mode feedback is 
active only in massive objects and at late times and it has 
no effect during the principal growth phase of most galaxies 
and AGN. 

The distribution of galaxy Ught in our 'best' model is 
shown in the bottom panel of Fig. for comparison with 
the mass distribution in the top panel. For both the vol- 
ume is a projected 330 x 280 x 15/i~^Mpc slice cut from 
the full 0.125 /i~^Gpc"^ simulation box. The plot of surface 
brightness is colour-coded by the luminosity-weighted mean 
B — V colour of the galaxies. On large scales light clearly fol- 
lows mass, but non-trivial biases become evident on smaller 
scales, especially in 'void' regions. The redder colour of 
galaxies in high density regions is also very clear. 



3.2 Gas infall and cooling 

We fo llow the standard paradigm set out bv lWhite fc FrenkI 
il99ll) as adapted for implementation on high resolu- 
tion N-body s i mulat ions by [SprinEcl ct al. (20o2) ^^nd 
iDe Lucia et al.l (|20o3) ■ This assumes that as each dark mat- 
ter halo collapses, its own 'fair share' of cosmic baryons 
collapse with it (but see Section 13.31 below) . Thus in our 
model the mass fraction in baryons associated with every 
halo is taken t o be ff, = 17%, cons istent with the first-year 
WMAP result dSoergel et alJl2003ll . Initially these baryons 
are in the form of diffuse gas with primordial composition, 
but later they include gas in several phases as well as stars 
and heavy elements. The fate of the infalling gas depends on 



redshift and on the depth of the halo potential ('Silk"l977|; 
iRees fc OstrTkeilllQ??!: lBinnevlll977t IWhite fc Reea .197^) . 
At late times and in massive systems the gas shocks to the 
virial temperature and is added to a quasi-static hot atmo- 
sphere that extends approximately to the virial radius of the 
dark halo. Gas from the central region of this atmosphere 
may accrete onto a central object through a cooling flow. 
At early times and in lower mass systems the infalling gas 
still shocks to the virial temperature but its post-shock cool- 
ing time is sufficiently short that a quasi-static atmosphere 
cannot form. Rather the shock occurs at much smaller ra- 
dius and the shocked gas cools rapidly and settles onto a 
central object, which we assume to be a cold gas disk. This 
may in turn be subject to gravitational instability, leading 
to episodes of star formation. 

More specifically, the cooling time of a gas is conven- 
tionally taken as the ratio of its specific thermal energy to 
the cooling rate per unit volume, 



3 fiTUp kT 
^cool 77 " 



2p<,(r)A(r, Z) 



(3) 



Here fnUp is the mean particle mass, k is the Boltzmann 
constant, pg{r) is the hot gas density, and A{T,Z) is the 
cooling function. The latter depends both on the metallicity 
Z and the temperature of the gas. In our models we assume 
the post-shock temperature of the infalling gas to be the 
virial temperature of the halo, T = 35.9 ('K/ir/kms~^)^K. 
When needed, we assume that the hot gas within a static 
atmosphere has a simple 'isothermal' distribution. 



Pair) 



(4) 



where mhot is the total hot gas mass associated with the 
halo and is assumed to extend to its virial radius 7?vir. 

To estimate an instantaneous cooling rate onto the cen- 
tral object of a halo, given its current hot gas content, we 
define the cooling radius, rcooi , as the radius at which the lo- 
cal cooling time (assuming the struct ure of Eq.|l]l is equal to 
a suitably defined age for the halo. IWhite fc FrenkI 

119911) 

took this ag e to be the Hubble time fa. while Icol^t^ll 
il994l l200Cl) used the time s cale over which the ma i n pro - 
genitor last doubled its mass. ISomerville fc Primackl lll999h 
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Figure 2. The ratio of the cooling radius to virial radius for a random selection of viriahsed systems at 2 = plotted against their dark 
matter virial mass. Systems identified to be in the 'rapid cooling' regime are shown in the left panel, while those that ha ve formed a 
static hot halo are shown on the right (Section l3.2l . A sharp transition between the two regimes is seen close to that found bv lKeres et alj 
1200^^) ■ marked by the solid vertical line. 



argued that the time interval since the last major merger 
is more appropriate since such mergers redistribute hot gas 
within the halo. Here we follow Springel et al. (200l3) and 
iDe Lucia et al.l ll2004) and define the cooling radius as the 
point where the local cooling time is equal to the halo dy- 
namical time, -Rvir/Vvir = 0.1 H{z)~^ . This is an order of 
magnitude smaller than tn and so results in substantially 
smaller cooling radii and cooli ng rates (typically by a factor 
of 3) than the assumption of IWhite fc FrenkI Jl99ll) . Our 
choice is justified by the tests of ^foshM^^^d] 12002) who 
verified explicitly that it results in good object-by-object 
agreement between the amount of gas predicted to con- 
dense in galaxy mass haloes and the amount which actually 
condensed in their high-resolution SPH simulations of the 
formation of a galaxy cluster and its environment. These 
tests assumed primordial abundances in the cooling func- 
tion. When we implement a chemical enrichment model con- 
sistent with the observed element abundances in intracluster 
gas (see Section r3.9t cooling rates in galaxy mass haloes are 
substantially enhanced and we find (as did IPe Lucia et all 
[2004) that a sm aller coefficient than used in the original 
Iwhite fc FrenkI cooling model is required to avoid excessive 
condensation of gas. 

Using the above definition, a cooling rate can now be 
determined through a simple contiimity equation, 

mcooi = 'ii^pg(rcooi)r'Loi'rcooi ■ (5) 

Despite its simplicity, this equation is a good approximation 
to the rate a t whic h gas is deposited at the centre in the 
iBertschingeil I^S^) similarity solution for a cooling fiow. 
Putting it all together we take the cooling rate within a 
halo containing a hot gas atmosphere to be 

mcool = 0.5 rrihot — ^ ■ (6) 

We assume this equation to be valid when rco oi < -Rvir. 
This is the criterion which Iwhite fc FrenkI jlQQlll proposed 
to separate the static hot halo regime from the rapid cooling 
regime (see below). 

In low mass haloes or at high redshifts the formal cool- 



ing radius lies outside the virial radius rcooi > Rvir- The 
post-shock gas then cools in less than one sound crossing 
time and cannot maintain the pressure needed to support an 
accretion shock at larg e radius. The high-resolution s pheri- 
cal infall simulations of lForcada-Miro fc Whit^ I^Q^) show 
that in this situation the accretion shock moves inwards, 
the post-shock temperature increases and the mass stored in 
the post-shock hot atmosphere decreases, because the post- 
shock gas rapidly cools onto the central object. In effect, all 
infalling material is accreted immediately onto the central 
disk. In this rapid cooling regime we therefore set the cooling 
rate onto the central object to be equal to the rate at which 
new diffuse gas is added to the halo. 



3.2.1 Rapid cooling or cold accretion? 

Although much simp lifie d, the above mode l was shown by 
lYoshida et"ai] ll2002li and lHellv et alJ J2003t) to give reason- 
ably accurate, object-by-object predictions for the cooling 
and accumulation of gas within the galaxies that formed 
in their Af-body-|-SPH simulations. These neglected star- 
formation and feedback effects in order to test the cooling 
model alone. They also assumed primordial abundances in 
the cooling function. In Fig. |5|we show the ratio rcooi/i?vir 
as a function of virial mass for haloes in the 'rapid cool- 
ing' regime (left panel) and in the 'static halo' regime (right 
panel) at z = for our 'best' galaxy formation model. The 
two regimes are distinguished by the dominant gas com- 
ponent in each halo: when the mass of hot halo gas ex- 
ceeds that of cold disk gas, we say the galaxy has formed 
a static halo, otherwise the system is taken to be in the 
rapid cooling phase. Many haloes classified as 'rapidly cool- 
ing' by this criterion have rcooi < Rvir, which would appar- 
ently indicate a static hot halo. This is misleading, how- 
ever, as systems where cooling is rapid deposit infalling 
gas onto the central galactic disk on a short timescale 
until they have a low-mass residual halo which satisfies 
'"cool ~ Rvir. This then persists until the next infall event. 
Averaging over several cycles of this behaviour, one finds 
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that the bulk of the infalling gas cools rapidly. This is why 
we choose to classify systems by their dominant gas com- 
ponent. Note also that a massive hot halo forms imme- 
diately once cooling be comes inefficient, just as in t he 1- 
D infall simulations of | |brcada-Mir6 fc Whit^ (^O^) and 
iBirnboim fc DekeJ 1)200.^^ . Our classification is thus quite 
robust. 

The transition between the 'rapid cooling' and 'static 
halo' regimes is remarkably well defined. At z = it oc- 
curs at a halo virial mass of 2-3 x 10^^ M©, and is approx- 
imately independent of redshift out to at least z = 6. This 
is clos e to the transition mass f ound for the same cosmol- 
ogy by IBirnboim fc Dekel| j|200j)_using spherically symmet- 
ric simulations, and bv iKere^etal] (p004) using fully 3-D 
simulations. This is, in fact, a coincidence since both sets of 
authors assume cooling functions with no metals, whereas 
our 'best' model includes heavy elements which substantially 
enhance cooling at the temperatures relevant for the transi- 
tion. Enrichment in the 'rapid cooling' regime arises from the 
mixing of infalling primordial gas with supernova-enriched 
gas ejected in a galactic wind (see Section r3.6ll . If we modify 
our model to assume a zero-metal cooling function, we find 
the transition mass to shift downwards by about a factor of 
2-3, res ulting in a lower cooling rate in comparison to t he 
work of IBirnboim fc Dekell (1200311 and lKeres etall J2004l . 

The reason for the more efficient (we would argue overly 
efficient) cooling appears to be different in those two studies. 
The s pherical infall simulations of £orcada-Mir6 fc White! 
1^93) showed good agreement with a transition at the 
p oint predicted us i ng th e original cooling radius definition 
of I White fc FrenkI (119911 ) rather than our revised definition 
which was checked explicitly against SPH simulations by 
FYoshida c t al. ( 20021). Spherical models thus predict more ef- 
ficient cooling than o ccurs in typical 3-D situations. This ex- 
plains, perhaps, whv lBirnboim fc DekeJ (|200^ find a higher 
tra nsition mass than we would predict for zero metallic- 
ity. lYoshida et al.l (1200 2 fl also showed that the density es- 
timati on in the implementation of SPH used bv lKeres et alJ 
l|20^ leads to overcooling of gas in galaxy mass objects as 
compared to their own entropy and energy con serving SPH 
schem e; the effect is large enough to explain whv lKeres et alJ 
i2004l) find a higher transition mass than we find for their 
assumed cooling function. 



Both IBirnboim fc Dekell and iKeres et al.l refer to the 
'rapid cooling' regime as 'cold infall'. This is, in fact, a mis- 
nomer. In this mode the accretion shock occurs closer to the 
central object, and so deeper in its potential well than when 
there is a static hot halo. As a result, the pre-shock veloc- 
ity of infalling gas is greater, resulting in a larger post-shock 
temperature. The two modes do not differ greatly in the tem- 
perature to which infalling gas is shocked, but rather in how 
long (compared to the system crossing time) the gas spends 
at the post-shock temperature before its infall energy is lost 
to radiation. Finally, we note that the existence and impor- 
tance of th ese two modes was the maj or in sight of the origi- 
nal wo rk of lsii^ (|l97i1) . lBinn'evl (119771) and lRees fc Ostrikeil 
il977tl and has been built into all modern theories for galaxy 
fo rmation. A detai l ed dis cussion can be found, for example, 
in lWhite fc FrenkI Jl99J) . 



3.3 Reionization 

Accretion and cooling in low mass haloes is required to 
be inefficient to explain why dwarf galaxies contain a rela- 
tively small fraction of all condensed baryons dWhite fc ReesI 
0^978). This inefficiency may in part result from photoionisa- 
tion heating of the intergalactic medium (IGM) which sup- 
presses the conce ntration of baryons in shallow potentials 
llEfstathioulll992l) . More recent models identify the possible 
low-redshift signature of such heating in the faint end of 
the galaxy luminosity function, most notably in the abun- 
dance of the dwarf satellite galaxies in the local group (e.g. 
iTuUy et al.ll2002l: iBenson et aljliooll . 

^~~ lGn^dir " il2000r ~showed that the effect of photoioniza- 
tion heating on the gas content of a halo of mass Mvir can 
be modelled by defining a characteristic mass scale, the so 
called filtering mass, Mp , below which the gas fraction ft is 
reduced relative to the universal value: 



/r°(^,MviO = 



/b 



• (7) 

(1 +0.26MF(2)/Mvir)3 ^ ' 

The filtering mass is a function of redshift and changes most 
significantly around the epoch of reionization. It was esti- 
mated by G nedin using high-r e solution SLH-P ^ M sim ula- 
tions (but see lHoeft et alJl2005l) . |Kravtsov et al.l f2004'l pro- 
vided an analytic model for these results which distinguishes 
three 'phases' in the evolution of the IGM: z> zo, the epoch 
before the first HII regions overlap; zq < z < Zr, the epoch 
when multiple HII regions overlap; z < Zr, the epoch when 
the medium is almost fully reionized. They find that choos- 
ing zo — S and Zr — J provides the best fit to the numerically 
determined filtering mass. We adopt these parameters and 
keep them fixed throughout our paper. This choice results in 
a filtering mass of 4 X 10^ M© at z = Zr, and 3 x 10^° Mp by 
the present day. See Appendix B of iKravtsov et all (|200J) 
for a full derivation and description of the analytic model. 



3.4 Black hole growth, AGN outflows, and 
cooling suppression 

There is a growing body of evidence that active galactic 
nuclei (AGN) are a critical piece in the galaxy formation 
puzzle. Our principal interest in this paper is their role in 
suppressing cooling flows, thereby modifying the luminosi- 
ties, colours, stellar masses and clustering of the galaxies 
that populate the bright end of the galaxy luminosity func- 
tion. To treat this problem, we first need a physical model 
for the growth of black holes within our galaxies. 



3.4-1 The 'quasar mode' 

In our model (which is based closely on that of 
Kauffmann & Haclmclt |200(J) supermassive black holes 
grow during galaxy mergers both by merging with each other 
and by accretion of cold disk gas. For simplicity, black hole 
coalescence is modelled as a direct sum of the progenitor 
masses and thus ignores gravitational wave losses (including 
such losses is found to have little effect on the properties of 
the final galaxy population). We assume that the gas mass 
accreted during a merger is proportional to the total cold 
gas mass present, but with an efficiency which is lower for 
smaller mass systems and for unequal mergers. Specifically, 
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AmsH.Q = 



/bH "Icold 



(8) 



1 + (280kms-VVvir)2 ' 

where we have changed the original parameterisation by tak- 
ing 

/bH = /bH (msat/mcontral) ■ (9) 

Here /bh ~ 0.03 is a constant and is c hosen to reproduce the 
observed local msH — n t buire relation llMagorrian et alJl998t 



Marconi fc Huntll2003l: 



Haring fc Rixll2004F " In contrast to 



Kauffmann fc Haehneltl 11200(1) we allow black hole accre- 
tion during both major and minor mergers although the ef- 
ficiency in the latter is lower because of the maat /n-contrai 
term. Thus, any merger-induced perturbation to the gas 
disk (which might come from a bar instability or a merger- 
induced starburst - see Section [3.711 can drive gas onto the 
central black hole. In this way, minor merger growth of the 
black hole parallels minor merger growth of the bulge. The 
fractional contribution of minor mergers to both is typically 
quite small, so that accretion driven by major mergers is the 
dominant mode of black hole growth in our model. We refer 
to this as the 'quasar mode'. [Note that a more schematic 
treatment of black hole gro wth would suffice for the pur- 
poses of this paper, but in ISpringel et all ^l2005b^ and in 
future work we wish to examine the build-up of the black 
hole population within galaxies in considerably more detail.] 
There is substantial evidence for strong hydrodynamic 
and radiative fee d back from optical/UV and X-ray AGN 
jArav et alJ l200ll: Ide Kool et alj lioOlt iReeves et all 120031: 
ICrenshaw et al.l 120031) .' We have not yet explicitly incor- 
porated such feedback in our modelling, and it may well 



turn out to 



for 



lie, the recent 



simulations of 


Di Matteo et alJl2005l: ISoringel et al.ll2005al: 


FlIoDkins et alj 


20051). We assume 'quasar mode' accretion 



to be closely associated with starbursts, so this feedback 
channel may be partially represented in our models by an 
enhanced effective feedback efficiency associated with star 
formation and supernovae (see Section [3.61 . 
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Figure 3. The black hole accretion rate density, jtibhi a 
function of redshift for both the 'quasar' and the 'radio' modes 
discussed in Section 13.41 This figure shows that the growth of 
black holes is dominated by the 'quasar mode' at high redshift and 
falls off sharply at 2 ^ 2. In contrast, the 'radio mode' becomes 
important at low redshifts where it suppresses cooling flows, but 
is not a significant contributor to the overall black hole mass 
budget. 



is typically orders-of-magnitude below the Eddington limit. 
In Section|5]we discuss physical accretion models which may 
reproduce this phenomenology. 

We assume that 'radio mode' feedback injects sufficient 
energy into the surrounding medium to reduce or even stop 
the cooling flow described in Section 13.21 We take the me- 
chanical heating generated by the black hole accretion of 
Eq.Hnito be 



JBH = 7?fnBH C 



(11) 



3.4-2 The 'radio mode' 

In our model, low energy 'radio' activity is the result of 
hot gas accretion onto a central supermassive black hole 
once a static hot halo has formed around the black hole's 
host galaxy. We assume this accretion to be continual and 
quiescent and to be described by a simple phenomenological 
model: 

msH \ ( ./"hot \ / Vvir 



JTlBH.R = KAGN 



( 108 Mo 



(tr)( 



200 kms- 



(10) 



where ttibh is the black hole mass, /hot is the fraction of 

1/2 

the total halo mass in the form of hot gas, K/ir oc TJ^^ is 
the virial velocity of the halo, and kagn is a free parameter 
with units of Moyr^^ with which we control the efficiency 
of accretion. We find below that kagn = 6 x lO~®M0yr~^ 
accurately reproduces the turnover at the bright end of the 
galaxy luminosity function. Note that /hotKfir^H is propor- 
tional to the total mass of hot gas, so that our formula is 
simply the product of the hot gas and black hole masses mul- 
tiplied by a constant efficiency and divided by the Hubble 
time tn. In fact, we find /hot to be approximately constant 
for Vvir ~ 150kms~^, so the dependence of jtibh.r on this 
quantity has little effect. The accretion rate given by Ea. llOl 



where r\ = 0.1 is the standard efficiency with which mass 
is assumed to produce energy near the event horizon, and c 
is the speed of light. This injection of energy compensates 
in part for the cooling, giving rise to a modified infall rate 
(Eq.il of 



'^cool " ^^cool 



2 vir 



(12) 



For consistency we never allow ra^^^^ to fall below zero. It is 



worth noting that rhcooi oc A(Vvir) 



1/2 



'1/2 



(Eq.EJ 



and riihoat = 2Lbii/VX « mBR /hot Kir (Eq.^, so that 



W2heat 
?^cool 



,1/2 

rriBH 



/ho^t A(Kir)V2Kir 



(13) 



(These scalings are exact for an EdS universe; we have omit- 
ted weak coefficient variations in other cosmologies.) Thus 
in our model the effectiveness of radio AGN in suppressing 
cooling flows is greatest at late times and for large values 
of black hole mass. This turns out to be the qualitative be- 
haviour needed for the suppression of cooling flows to repro- 
duce successfully the luminosities, colours and clustering of 
low redshift bright galaxies. 
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Figure 4. The black hole-bulge mass relation for model galaxies 
at the present day. The local observational result of lHaring fc RiJ 
JSooJ is given by the solid line, where the dashed box shows the 
approximate range over which their fit was obtained. 



3.4-3 The growth of supermassive black holes 

Fig. 121 shows the evolution of the mean black hole accretion 
rate per unit volume averaged over the entire Millennium 
Simulation. We separate the accretion into the 'quasar' and 
'radio' modes described above (solid and dashed lines re- 
spectively). Black hole mass growth in our model is dom- 
inated by the merger-driven 'quasar mode', which is most 
efflcient at redshifts of two to four, dropping by a factor 
of five by redshift zero. This behaviour has similar form to 
but is weaker than the observed evolution with redshift o f 
the bright quasar population (a^jH^rtwig^^^chad^ 1 99d|) • 
(See also the discussion in Kauf fmann fc Haehn elt 2000). In 
contrast, our 'radio mode' is significant only at late times, as 
expected from the scaling discussed above, and for the high 
feedback efficiency assumed in Eg. llll it contributes only 5% 
of the final black hole mass density. We will show, however, 
that the outflows generated by this accretion can have a 
major impact on the final galaxy properties. Finally, inte- 
grating the accretion rate density over time gives a present 
day black hole mass density of 3 x Mp,Mx)C~^ , consistent 
with recent observational estimates jYu fc Tremainejl2002t 
^le rlon i 2004). 

The relationship between black hole mass and bulge 
mass is plotted in Fig. |1] for the local galaxy population in 
our 'best' model. In this figure, t he solid line shows th e best 
fit to the observations given bv lHarin~ fc Rixl J2004 for a 
sample of 30 nearby galaxies with well measured bulge and 
black hole masses. Their results only probe masses over the 
range bounded by the dashed lines. Our model galaxies pro- 
duces a good match to these observations with comparable 
scatter in the observed range (see their Fig. 2). 



3.5 Star formation 

We use a simple model for star formation similar to those 
adopted by earlier authors. We assume that all star for- 
mation occurs in cold disk gas, either quiescently or in a 
burst (see Sectio n 13.71 . Based on the observational work of 
iKennicutJ ilQQST) . we adopt a threshold surface density for 
the cold gas below which no stars form, but abo ve which gas 
starts to collapse and form stars. According to iKauffmannI 
il99(f), this critical surface density at a distance 7? from the 
galaxy centre may be approximated by 

We convert this critical surface density into a critical mass 
by assuming the cold gas mass to be evenly distributed over 
the disk. The resulting critical cold gas mass is: 

— 3.8X10^ (^^^)(^)m,, (15) 

where we assum e the disk scale length to be Ts = (A/\/2)i?vir 
||Mo et al. ^l998l. and set the outer disk radiu s to rdisk = 3rj, , 
based on the properties of the Milky Way (Ivan den Berghl 
|2Q00'') • Here A is the spin parameter of th e dark halo in which 
the galaxy resides ijBuUock et ahlEoOlT l. measured directly 
from the simulation at each time step. When the mass of 
cold gas in a galaxy is greater than this critical value we 
assume the star formation rate to be 

rh* = asF (mcoid — merit) / tdyn,disk , (16) 

where the efficiency asF is typically set so that 5 to 15 per- 
cent of the gas is converted into stars in a disk dynamical 
time tdyn.disk, which we define to be rdisk/'K/ir. This star 
formation model produces a global star formation history 
consistent with the observed star formation density of the 
universe out to at least z = 2, as shown in Fig. |^ (note that 
this figure also includes star formation through starbursts - 
see Section f3.7ll . 

When implemented in our model, Eq. 1161 leads to 
episodic star formation that self-regulates so as to main- 
tain the critical surface density of Eq. 1141 This naturally 
reproduces the observed spiral galaxy gas fractions with- 
out the need for additional parameterisation, as we demon- 
strate in the top panel of Fig. |S| using model Sb/c galaxies 
identified as objects with bulge-to-total luminosity: 1.5 ^ 
Mi, bulge — Afi, total ^ 2.5 (bulge formation is described in 
Section 13.711 . 

3.6 Supernova feedback 

As star formation proceeds, newly formed massive stars 
rapidly complete their evolution and end their life as super- 
novae. These events inject gas, metals and energy into the 
surrounding medium, reheating cold disk gas and possibly 
ejecting gas even from th e surrounding halo. 

The observations of iMartinI lll999l) suggest modelling 
the amount of cold gas reheated by supernovae as 

Arrircheatcd = ediskAm* , (17) 

where Am* is the mass of stars formed over some finite time 
interval and tdisk is a parameter which we fix at edisk = 3.5 
based on the observational data. The energy released in this 
interval can be approximated by 
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Figure 5. The star formation rate density of the universe as a 
function of redshift. The symbols show a compilation of observa- 
tional results taken from Fig. 12 of Soringcl & Hcrnquist 1 20033)- 
The solid line shows our 'best' model, which predicts that galax- 
ies form much of their mass relatively early. The dashed line (and 
right axis) indicate the increase in stellar mass with redshift. Ap- 
proximately 50% of all stars form by z = 2. 



AEs 



0.5 Ehaio Am*Vs] 



(18) 



where 0.5 Vg^ is the mean energy in supernova ejecta per 
unit mass of stars formed, and ehaio parametrises the effi- 
ciency with which this energy is able to reheat disk gas. 
Based on a standard initial stellar mass function and stan- 
dard supernova theory we take Vsn ~ 630kms~^. In ad- 
dition, for our 'best' model we adopt ehaio = 0.35. If the 
reheated gas were added to the hot halo without changing 
its specific energy, its total thermal energy would change by 

^2 



AEkot = 0.5 Am„ 



(19) 



Thus the excess energy in the hot halo after reheating is just 
ASexccss = A_BsN - A^hot- When ABcxcoss <0 the energy 
transferred with the reheated gas is insufficient to eject any 
gas out of the halo and we assume all hot gas remains as- 
sociated with the halo. When excess energy is present, i.e. 
Ai?oxcess >0, we assume that some of the hot gas is ejected 
from the halo into an external 'reservoir'. Specifically, we 
take 



A?71ejected 



rrihot 



(ehaio^-ediBk)Am. , (20) 



-Ehot 

where Shot = 0.5 mhotK^ir is the total thermal energy of the 
hot gas, and we set Amejoctcd = when this equation gives 
negative values (implying A_Ecxccss <0 as discussed above). 
This is similar to the traditional semi-analytic feedback 
recipe, Amejoctcd oc Arrit/V^i^, but with a few additions. 
For small K-ir the entire hot halo can be ejected and then 
Amejectcd must saturate at Arrireheated- Conversely, no hot 
gas can be ejected from the halo for V^^^ > ehaioV's^N/^disk, 
i.e. for halo circular velocities exceeding about 200kms~^ 
for our favoured parameters. 

Ejected gas leaves the galaxy and its current halo in a 
wind or 'super-wind', but it need not be lost permanently. 
As the dark halo grows, some of the surrounding ejecta may 
fall b ack in and be re i ncorpo rate d into the cooling cycle . We 
foUow lSprineel et al.l J2001all and lDe Lucia et all J2004) and 
model this by assuming 
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Figure 6. Selected results for Sb/c galaxies (identified by bulge- 
to-total luminosity, see Section l3.5i for our best model, (top) Gas 
fractions as a function of B-band magnitude. The solid line is a 
representation of the mean behaviour in the (incomplete) sample 



of lGarnetj j2002). (middle) The TuUy-Fisher relation, where the 
disk circular velocity, Vc, is approximated by Vvir for the dark 
halo. The solid line with surrounding dashed lines represents the 
mean result and scatter found by jGiovanelli et al. 11997). (bot- 
tom) Cold gas metallicity as a fu nction of total stellar m ass. The 
solid line represents the result of lTremonti et aP J2004l) . 
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^ejected — O'cj ^ejected / ^dyn j (*^-^) 

where 7cj controls the amount of reincorporation per dy- 
namical time; typically we take 7oj = 0.3 to 1.0. Such values 
imply that all the ejected gas will return to the hot halo in 
a few halo dynamical times. 

The prescriptions given in this section are simple, as 
well as physically and energetically plausible, but they have 
little detailed justification either from observation or from 
numerical simulation. They allow us to track in a consis- 
tent way the exchange of each halo's baryons between our 
four phases (stars, cold disk gas, hot halo gas, ejecta), but 
should be regarded as a rough attempt to model the complex 
astrophysics of feedback which will surely need significant 
modification as the observational phenomenology of these 
processes is explored in more depth. 

Supernova feedback and star formation act together to 
regulate the stellar growth of a galaxy. This is especially im- 
portant for L < L* galaxies, where feedback can eject most 
of the baryons from the system, reducing the supply of star- 
forming material for time periods much longer than the cool- 
ing/supernova heating cycle. In the middle panel of Fig. |H] 
we plot the TuUy-Fisher relation for model Sb/c galaxies 
(see Section 13.511 . The TuUy-Fisher relation is strongly in- 
fluenced by the link between star formation and supernova 
heating. The circular velocity of a galactic disk is (to first 
order) proportional to the virial velocity of the host dark 
matter halo and thus to its escape velocity. In our model 
(and most others) this is closely related to the ability of the 
galaxy to blow a wind. The luminosity of a galaxy is deter- 
mined by its ability to turn its associated baryons into stars. 
The overall efficiency of this process in the face of supernova 
and AGN feedback sets the amplitude of the TuUy-Fisher re- 
lation, while the way in which the efficiency varies between 
systems of different circular velocity has a strong influence 
on the slope. 

To predict a TuUy-Fisher relation for our model we need 
to assign a maximum rotation velocity to each of our galaxy 
disks. For central galaxies we simply take this velocity to be 
the Vvir of the surrounding halo, while for satellite galax- 
ies we take it to be the Vvir of the surrounding halo at the 
last time the galaxy was a central galaxy. This is a crude 
approximation, and for realistic halo structures it is likely 
to be an underestimate both because of the concentration 
of the dark matter distribution and because of the grav- 
itational effects of the baryons (e.g. iMo et ahl 1199^1 . Ob- 
taining a good simultaneous flt to observed luminosity func- 
tions and TuUy-Fisher relations remains a diffic ult problem 
withi n the ACDM paradigm (see, for example ICole et alj 
|200[I) . Our unrealistic assumption for the disk rotation ve- 
locity actually produces quite a good fit to the observational 
data of lGiovanelli et al.l lll997|l . demonstrating that our sim- 
ple star formation and feedback recipes can adequately rep- 
resent the growth of stellar mass across a wide range of 
scales. We find clear deviations from power law behaviour 
for log VK ^ 2.3 (approximately Vc Si lOOkms"^), where the 
efficiency of removing gas from low mass systems combines 
with our threshold for the onset of star formation to reduce 
the number of stars that can form. The resulting downward 
be nd is qualitatively simil ar to that pointed out in real data 
bv lMcGaugh et all i2000l) . These authors show that includ- 
ing the gaseous component to construct a 'baryonic' TuUy- 



Fisher relation brings the observed points much closer to a 
power-law, and the same is true in the model we present 
here. Limiting star formation in galaxies that inhabit shal- 
low potentials has a strong effect on the faint-end of the 
galaxy luminosity function, as will be seen in Section. 14.21 

3.7 Galaxy morphology, merging and starbursts 

In the model we discuss here, the morphology of a galaxy 
is assumed to depend only on its bulge-to-total luminosity 
ratio, which in turn is determined by three distinct physical 
processes: disk growth by accretion, disk buckling to produce 
bulges, and bulge formation through mergers. We treat disk 
instabiliti es using the simple analytic stability criterion of 
IMo et al] (1998); the stellar disk of a galaxy becomes unsta- 
ble when the following inequality is met. 



(GmD/rD)i/2 " ' ^ ^ 

where we again approximate the rotation velocity of the disk 
Vc by Vvir. For each galaxy at each time-step we evaluate the 
left-hand side of Eq. 1221 and if it is smaller than unity we 
transfer enough stellar mass from disk to bulge (at fixed ru) 
to restore stability. 

Galaxy mergers shape the evolution of galaxies, affect- 
ing both their morphology and (through induced starbursts) 
their star formation history. Mergers can occur in our model 
between the central galaxy of a dark halo or subhalo and a 
satellite galaxy which has lost its own dark subhalo. Sub- 
structure is followed in the Millennium Run down to a 
20 particle limit, which means that the orbit of a satellite 
galaxy within a larger halo is followed explicitly until its sub- 
halo mass drops below 1.7 x 10^*^ Mq. After this point, 
the satellite's position and velocity are represented by those 
of the most bound particle of the subhalo at the last time it 
was identified. At the same time, however, we start a merger 
'clock' and estimate a merging ti me for the galaxy using the 
dynamical friction formula of Binnev fc Tremaind (^83), 

tMction = 1-17 /"""^T,. ■ (23) 

Grrisat In A 

This formula is valid for a satellite of mass rriaat orbiting 
in an isothermal potential of circular velocity Vvir at ra- 
dius Tsat. We take msat and rsat to be the values mea- 
sured for the galaxy at the last time its subhalo could 
be identified. The Coulomb logarithm is approximated by 
InA = ln(l -I- Mvir/jTisat). The satellite is then merged with 
the central galaxy a time tfriction after its own subhalo was 
last identified. If the main halo merges with a larger system 
before this occurs, a new value for ttriction is calculated and 
the merger clock is restarted. 

The outcome of the merger will depend on the baryonic 
mass ratio of the two progenitors . When one dominates the 
process, i.e. a small satellite merging with a larger central 
galaxy, the stars of the satellite are added to the bulge of the 
central galaxy and a minor merger starburst (see below) will 
result. The cold gas of the satellite is added to the disk of 
the central galaxy along with any stars that formed during 
the burst. Such an event is called a minor merger. 

If, on the other hand, the masses of the progenitors are 
comparable a major merger will result. Under these circum- 
stances the starburst is more significant, with the merger 
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Figure 7. The mean condensation rate, (riieool) ^.s a function of halo virial velocity V^ij at redshifts of 6, 3, 1, and 0. Solid and dashed 
lines in each panel represent the condensation rate with and without 'radio mode' feedback respectively, while the vertical dotted lines 
show the transition between the rapid cooling and static hot halo regimes, as discussed in Section i3.2l This figure demonstrates that 
cooling flow suppression is most efficient in our model for haloes with V^ir > 150kms~^ and at 2 ^ 3. 



destroying the disks of both galaxies to form a spheroid in 
which all stars are placed. The dividing line between a ma- 
jor and minor merger is given by the parameter Tmerger: 
when the mass ratio of the merging progenitors is larger 
than Tincrgcr a major mer ger results, otherwise t he event is 
a minor merger. Following ISpringel et al.l i2001al) we choose 
Tmcrgcr = 0.3 and keep this fixed throughout. 

Our starburst implementation is based on the 'coUi- 
sional starburst' model of ^Spmorvillo ot al. (2001). In this 
model, a fraction eburst of the combined cold gas from the 
two galaxies is turned into stars as a result of the merger: 

eburst = /3burat(msat/mcentral , (24) 

where the two parameters are taken as aburst = 0.7 and 
/3burst = 0.56 . This mod e l prov ides a go od fit to the numeri- 
cal results o f ICox et alJ i2004) and also lMihos fc Hernauisd 
for merger mass ratios ranging from 1:10 to 1:1. 



3.8 Spectroscopic evolution and dust 

The photometric properties of our galaxies are cal- 
culated using stel l ar po pulation synthesis models from 
iBruzual fc Charlol] l|20Q3'l. Our implementation is fully de- 
scribed "in TD^Lu^i^^^ l. (2004) and we refer the reader 
there (and to references therein) for further details. 

To include the effects of dust when calculating galaxy 
lu minosities we apply the s imple 'plane-parallel slab' model 
of iKauffmann et al.l lll999l) . This model is clearly oversim- 
plified, but it permits us to make a reasonable first-order 
correction for dust extinction in actively star-forming galax- 
ies. For the deta i ls of this model we refer the reader to 
IKauffmann et akl (Il999l) and to references therein. 

3.9 Metal enrichment 

Our treatment of metal enrichment is ess entially identical 
to that described in lDe Lucia et aP (|20o3). In this model a 
yield Y of heavy elements is returned for each solar mass of 




Figure 8. Galaxy luminosity functions in the K (left) and bj (right) photometric bands, plotted with and without 'radio mode' feedback 
(solid and long dashed lines respectively — see Section l3.4l . Symbols indicate observational results as listed in each panel. As can be seen, 
the inclusion of AGN heating produces a good fit to the data in both colours. Without this heating source our model overpredicts the 
luminosities of massive galaxies by about two magnitudes and fails to reproduce the sharp bright end cut-offs in the observed luminosity 
functions. 



stars formed. These metals are produced primarily in the su- 
pernovae which terminate the evolution of short-lived, mas- 
sive stars. In our model we deposit them directly into the 
cold gas in the disk of the galaxy. (An alternative would 
clearly be to add some fraction of the metals directly to 
the hot halo. Limited experiments suggest that this makes 
little difference to our main results.) We also assume that 
a fraction R of the mass of newly formed stars is recycled 
immediately into the cold gas in the disk, the so called 'in - 
stantaneous recycling approximation' (see lCole et aLlEoOCtl . 
For full details on metal enrichment and exchange processes 
in our model seejOc Lucia et al. (2004). In the bottom panel 
of Fig. El we show the metallicity of cold disk gas for model 
Sb/c galaxies (selected, as before, by bulge-to-total luminos- 
ity, as described in Section FS.St as a functio n of total stellar 
mass. For comparison, we show the result of iTre monti et all 
l|2004l for mean HII region abundances in SDSS galaxies. 



4 RESULTS 

In this section we examine how the suppression of cooling 
flows in massive systems affects galaxy properties. As we will 
show, the effects are only important for high mass galaxies. 
Throughout our analysis we use the galaxy formation model 
outlined in the previous sections with the parameter choices 
of Table except where explicitly noted. 



4.1 The suppression of cooling flows 

We begin with Fig. [3 which shows how our 'radio mode' 
heating model modifies gas condensation. We compare mean 
condensation rates with and without the central AGN heat- 
ing source as a function of halo virial velocity (solid and 
dashed lines respectively). Recall that virial velocity pro- 
vides a measure of the equilibrium temperature of the sys- 
tem through Tvir oc V^i,., as indicated by the scale on the top 
axis. The four panels show the behaviour at four redshifts 
between six and the present day. The vertical dotted line in 
each panel marks haloes for which rcooi ~ Rvir, the transi- 
tion between systems that form static hot haloes and those 
where infalling gas cools rapidly onto the central galaxy disk 
(see section Isl and Fig.|2j. This transition moves to haloes 
of lower temperature with time, suggesting a 'down-sizing' of 
the characteristic mass of actively star-forming galaxies. At 
lower Vvii gas continues to cool rapidly, while at higher Vvir 
new fuel for star formation must come from cooling flows 
which are affected by 'radio mode' heating. 

The effect of 'radio mode' feedback is clearly substan- 
tial. Suppression of condensation becomes increasingly effec- 
tive with increasing virial temperature and decreasing red- 
shift. The effects are large for haloes with Vvir ^ 150kms~^ 
(Tvir lO^K) at 2^3. Condensation stops almost com- 
pletely between 2 = 1 and the present in haloes with 
Kir > 300kms"^ (Tvir > 3 x lO'^K). Such systems corre- 
spond to the haloes of groups and clusters which are typ- 
ically observed to host massive elliptical or cD galaxies at 
their centres. Our scheme thus produces results which are 
qualitatively similar to the ad hoc suppression of cooling 
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Figure 9. The B — V colours of model galaxies plotted as a func- 
tion of their stellar mass with (top) and without (bottom) 'radio 
mode' feedback (see Section l3.4i . A clear bimodality in colour is 
seen in both panels, but without a heating source the most mas- 
sive galaxies are blue rather than red. Only when heating is in- 
cluded are massive galaxies as red as observed. Triangles (red) and 
circles (blue) correspond to early and late morphological types re- 
spectively, as determined by bulge-to-total luminosity ratio (see 
Section l4.2i . The thick dashed lines mark the resolution limit to 
which morphology can be reliably determined in the Millennium 
Run. 

flows ass umed in previo us models of galaxy formation. For 
example. lKauffmann et al. ( 1999) switched off gas condensa- 
tion in all haloes with K-ir > 350kms~^, while lHatton et all 
l)2fl03h stopped condensation when the bulge mass exceeded 
a critical threshold. 

4.2 Galaxy properties with and without AGN 
heating 

The suppression of cooling flows in our model has a dramatic 
effect on the bright end of the galaxy luminosity function. 
In Fig. |H| we present K- and bj-band luminosity functions 
(left and right panels respectively) with and without 'ra- 
dio mode' feedback (solid and dashed lines respectively). 
The luminosities of bright galaxies are reduced by up to 
two magnitudes when the feedback is switched on, and this 
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Figure 10. Mean stellar ages of galaxies as a function of stellar 
mass for models with and without 'radio mode' feedback (solid 
and dashed lines respectively). Error bars show the rms scatter 
around the mean for each mass bin. The suppression of cooling 
flows raises the mean age of high-mass galaxies to large values, 
corresponding to high formation redshifts. 



induces a relatively sharp break in the luminosity function 
which matches the observations wel l. We demonstrat e this 
by overplott i ng K- band data from ICole et all (1200 J) and 
Huang et alJ i20 0j) in the left panel, and bj-band data from 
Norberg et al.l (2002) in the right panel. In both band-passes 
the model is quite close to the data over the full observed 
range. We comment on some of the remaining discrepancies 
below. 

Our feedback model also has a significant effect on 
bright galaxy colours, as we show in Fig.|^ Here we plot the 
B — V colour distribution as a function of stellar mass, with 
and without the central heating source (top and bottom pan- 
els respectively). In both panels we have colour-coded the 
galaxy population by morphology as estimated from bulge- 
to-total luminosity ratio (split at I/buigc/itotai = 0.4). Our 
morphological resolution limit is marked by the dashed line 
at a stellar mass of ~ 4 x W^Mq; this corresponds approx- 
imately to a halo of 100 particles in the Millennium Run. 
Recall that a galaxy's morphology depends both on its past 
merging history and on the stability of its stellar disk in our 
model. Both mergers and disk instabilities contribute stars 
to the spheroid, as described in Section TS.?! The build-up of 
haloes containing fewer than 100 particles is not followed in 
enough detail to model these processes robustly. 

A number of important features can be seen in Fig. 1^1 
Of note is the bimodal distribution in galaxy colours, with 
a well-defined red sequence of appropriate slope separated 
cleanly from a broader 'blue cloud'. It is significant that 
the red sequence is composed predominantly of early-type 
galaxies, while the blue cloud is comprised mostly of disk- 
dominated systems. This aspect of our model suggests that 
that the physical processes which determine morphology 
(i.e. merging, disk instability) are closely related to those 
which control star formation history (i.e. gas supply) and 
thus determine galaxy colour. The red and blue sequences 
both display a strong metallicity gradient from low to high 
mass (c.f. Fig. ISJ, and it is this which induces a 'slope' in 




Figure 11. The bj-band galaxy luminosity function split by colour at B — V = 0.8 (Fig.|5J into blue (left panel) and red (right panel) sub- 
populations (solid lines ) . The dotted line s in ea ch panel repeat the opposite colour luminosity function for reference. Symbols indicate the 
observational results of lMadgwick et al.l J2002) for early and late-type 2dFGRS galaxies, split according to spectral type. Although our 
model split by colour captures the broad behaviour of the observed type-dependent luminosity functions, there are important differences 
which we discuss in Section 14.21 



the colour-magn itude relations w hich agrees well with ob- 
servation (e.g. B aldrv et"ani2004^ . 

By comparing the upper and lower panels in Fig. |U]we 
can see how 'radio mode' feedback modifies the luminosities, 
colours and morphologies of high mass galaxies. Not surpris- 
ingly, the brightest and most massive galaxies are also the 
reddest and are ellipticals when cooling flows are suppressed, 
whereas they are brighter, more massive, much bluer and 
typically have disks if cooling flows continue to supply new 
material for star formation. AGN heating cuts off the gas 
supply to the disk from the surrounding hot halo, truncating 
star formation and allowing the existing stellar population 
to redden. However, these massive red galaxies do continue 
to grow through merging. This mechanism allows the dom- 
inant cluster galaxies to gain a factor 2 or 3 in mass with- 
out significa nt star formation, in apparent agreement with 
observation ijAraeon-Salamanca et al. *1998). This late-stage 
(i.e. z ^ 1) hierarchical growth moves objects to higher mass 
without changing their colours. 

It is also interesting to examine the effect of AGN heat- 
ing on the stellar ages of galaxies. In Fig. 1101 the solid and 
dashed lines show mean stellar age as a function of stel- 
lar mass for models with and without 'radio mode' feed- 
back, while error bars indicate the rms scatter around the 
mean. Substantial differences are seen for galaxies with 
Afstoiiar ~ 10^^ Mq: the mean age of the most massive galax- 
ies approaching 12 Gyr when cooling flows are suppressed 
but remaining around 8 Gyr when feedback is switched off. 
Such young ages are clearly inconsistent with the old stel- 



lar populations observed in the majority of massive cluster 
ellipticals. 

The colour bimodality in Fig.|U]is so pronounced that it 
is natural to divide our model galaxies into red and blue pop- 
ulations and to study their properties separately. We do this 
by splitting at B — V = 0.8, an arbitrary but natural choice. 
Fig. llll shows separate bj luminosity functions for the result- 
ing populati ons. For comparison we overplot observational 
results from iMadgwick" et al.l (|2003) for 2dFGRS galaxies 
split by spectral type. Their luminosity functi ons are es- 
sentially identical to those of ICole et al.l ||^Q05) , who split 
the 2dFGRS catalogue by bj-rp colour. It thus can serve to 
indicate the observational expectations for populations of 
different colour. The broad behaviour of the red and blue 
populations is similar in the model and in the 2dFGRS. The 
faint-end of the luminosity function is dominated by late- 
types, whereas the bright end has an excess of early-types. 
The two populations have equal abundance about ha lf to 
one magnitude brighter than M^^ jNorberg et alj|2002t) . 

Fig. Illl also shows some substantial differences between 
model and observations. The red and blue populations differ 
more in the real data than they do in the model. There is a 
tail of very bright blue galaxies in the model, which turn out 
to be objects undergoing strong, merger-induced starbursts. 
These correspond in abundance, star formation rate and evo- 
lutionary state to the observed population of Ultraluminous 
Infrared Galaxies (ULIRG's) with the important difference 
that almost all the luminosity from young stars in the real 
systems is a bsorbed by dust and re-e mitted in the mid- to 
far-infrared dSanders fc Mirabellll99d) . Clearly we need bet- 
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ter dust modelling than our simple 'slab' model (Section lS.SH 
in order to reproduce the properties of such systems ade- 
quately. If we suppress starbursts in bright galaxy mergers 
we find that the blue tail disappears and the observed be- 
haviour is recovered. A second and substantial discrepancy 
is the apparent overproduction of faint red galaxies in our 
model, as compared to the 2dF measurem ents (however see 
iPopesso et 811120051 : iGonzalez era!] 12005*). Further work is 
clearly needed to understand the extent and significance of 
this difference. 



5 PHYSICAL MODELS OF AGN FEEDBACK 

Our phenomenological model for 'radio mode' feedback (Sec- 
tion 1^^^ is not grounded in any specific model for hot gas 
accretion onto a black hole or for the subsequent generation 
and thermalization of energy through radio outflows. Rather 
it is based on the observed properties of cooling flows and 
their central radio sources, and on the need for a source of 
feedback which can suppress gas condensation onto massive 
galaxies without requiring the formation of new stars. We 
have so far focused on the effects of such feedback without 
discussing how it might be realised. In this section we present 
two physical models which suggest how accretion onto the 
central black hole may lead to activity in a way which could 
justify the parameter scalings we have adopted. 

5.1 Cold cloud accretion 

A simple picture f or cooling flow evolut ion, based on the sim- 
ilarity solution of iBertschingeJ (^8^ for an unperturbed 
halo in isolation, can be summarised as follows. Cooling 
flows develop in any halo where the cooling time of the bulk 
of the hot gas is longer than the age of the system so that a 
static hot halo can form. Such haloes usually have a strong 
central concentration and we approximate their structure by 
a singular isothermal sphere. The inner regions then have a 
local cooling time shorter than the age of the system, and 
the gas they contain radiates its binding energy and flows 
inwards. The flow region is bounded by the cooling radius 
f cool where the local cooling time is equal to the age of the 
system (see section . This radius increases with time as 
t^^^. As Bertschinger showed, the temperature of the gas 
increases by about 20% as it starts to flow inwards, and its 
density profile fiattens to pg oc r~^^^. Initially, the flow is 
subsonic and each gas shell sinks stably and isothermally 
in approximate hydrostatic equilibrium. As it sinks, how- 
ever, its inward velocity accelerates because its cooling time 
shrinks more rapidly than the sound travel time across it, 
and at the sonic radius, rsonic, the two become equal. At 
this point the shell goes into free fall, its temperature de- 
creases rapidl y and it may fragment as a result of thermal 
insta bility llCowie et al.lll980ll:lNulsenlll986l:IBalbus fc Sokeil 
^8^. The dominant component of the infalling gas is then 
in the form of cold clouds and is no longer self-coupled by 
hydrodynamic forces. Different clouds pursue independent 
orbits, some with pericentres perhaps orders of magnitude 
smaller than rsonic- If these lie within the zone of infiuence of 
the black hole, tbh = GmBa/V^i^, we assume that some of 
the cold gas becomes available for fuelling the radio source; 
otherwise we assume it to be added to the cold gas disk. 



The parameter scalings implied by this picture can be 
estimated as follows. The sound travel time across a shell 
at the cooling radius is shorter than the cooling time by a 
factor ~ rcooi/-Rvir. At smaller radii the ratio of cooling time 
to sound travel time decreases as r^^'^ so that rsonic/rcooi ~ 
(rcooiZ-Rvir)^ implying rsonic ~ ''cooiZ-Rvir- If we adopt ren > 
10"* rsonic as the condition for effective fuelling of the radio 
source, we obtain 

mBH > 10"'' Afvir (rcool/-Rvir)^ (25) 

as the corresponding minimum black hole mass for frag- 
mented clouds to be captured. Under such conditions, only 
a small fraction (~ 0.01%) of the cooling flow mass need 
be accreted to halt the flow. The ratio in parentheses on 
the right-hand side of this equation scales approximately 
as rcooi/i?vir oc (mhot/Mvir)^''^iH^''^K7/' so ^he minimum 
black hole mass scales approximately as {m-tiot/Mvir)^^^t'^^^^ 
and is almost independent of Vvir- In our model, the growth 
of black holes through mergers and 'quasar mode' accretion 
produces a population where mass increases with time and 
with host halo mass. As a result, effective fuelling takes place 
primarily in the more massive haloes and at late times for 
this 'cold cloud' prescription. 

To test this particular model we switch off our stan- 
dard phenomenological treatment of 'radio mode' feedback 
(section 13.4.21 . assuming instead that feedback occurs only 
when Eo. 1251 is satisfied and that in this case it is sufficient 
to prevent further condensation of gas from the cooling flow. 
All other elements of our galaxy and black hole formation 
model are unchanged. The resulting cooling flow suppression 
is similar to that seen in Fig. |7] and all results presented 
in Section 121 and 0] are recovered. An illustration of this is 
given by Fig. 1121 where we compare the K-band luminos- 
ity function from this particular model (the dashed line) to 
the observational data (c.f. also Fig.jHJ. The model works so 
well, of course, because the numerical coefficient in Ea. l25l is 
uncertain and we have taken advantage of this to choose a 
value which puts the break in the luminosity function at the 
observed position. This adjustment plays the role of the effi- 
ciency parameter kagn in our standard analysis (see Ea. llUII . 

5.2 Bondi-Hoyle accretion 

Our second physical model differs from the first in assuming 
that accretion is not from the dominant, cold cloud compo- 
nent which forms within the sonic radius, but rather from a 
subdominant hot component which fills the space between 
these clouds. The clouds themselves are assumed to be lost 
to the star-forming disk. The dens ity profile of the r e sidua l 
hot component was estimated by iNulsen fc Fabiaiil ll200(f) 
from the condition that the cooling time of each radial shell 
should remain equal to the sound travel time across it as it 
flows inwards. This requires the density of the hot compo- 
nent to vary as 1/r within rsonic and thermal instabilities 
must continually convert material into condensed clouds in 
order to maintain this structure as the hot gas flows in. 

The rate at which hot gas is accreted onto the black 
hole can then be estimated from the Bondi-Hoyle formula 
(Bondi 1952; Edgar 200^: 

2 

niBondi = 2.57rG = — . (26) 
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Figure 12. The observed K-band galaxy luminosity function is 
compared with the results from models using the two physical pre- 
scriptions for 'radio mode' accretion discussed in Section |5] the 
Bondi-Hoyle accretion model (solid line) and the cold cloud ac- 
creti on model (dashed line ). Symbols indicate observational data 
from lCole et al.1 ^20nl^ and lHuang et al J 1200311 . Both models can 
produce a luminosity function which matches observation well. 

Here po is the (assumed uniform) density of hot gas around 
the black hole, and in all that follows we approximate the 
sound speed, Cg, by the virial velocity of the halo, Vvir- Of 
course, the density distribution of gas surrounding the black 
hole is not uniform so the question immediately arises as 
to what density we should choose. We follow a suggestion 
of E. Churazov and use th e value predicted by the 'maxi- 
mum cooling flow' model of iNulsen fc FabianI ll200dh at the 
Bondi radius rsondi = 2GA/bh/Cs = 2rBH, the conventional 
boundary of the sphere of influence of the black hole. We 
therefore equate the sound travel time across a shell at this 
radius to the local cooling time there: 

2rBon 



fimpkT 



Cs Ki 2pg(rBondi)A(r,Z) 

Solving for the density gives 



PO = /5g(rBondi 



kT Kir 



(27) 



(28) 
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Combining Eg . 1281 with 1261 provides us with the desired es- 
timate for the hot gas accretion rate onto the black hole: 



WlBondi 



kT 

G/imp — m-BH 



(29) 



Notice that this rate depends only on the black hole mass 
and on the virial temperature of the halo. It is independent 
both of time and of rrihot/A/vir, the hot gas fraction of the 
halo. It is valid as long as rBondi < »"sonic, which is always 
the case in our models. 

To investigate the effects of this model we replace the 
phenomenological 'radio mode' accretion rate of Eg. [Till with 
that given by Eq. 1291 Since the latter has no adjustable effi- 
ciency, we use the energy generation parameter rj of Eq. fTTI 
to control the effectiveness of cooling flow suppression. (This 
was not necessary before since rj always appeared in the 



product rj kagn , where kagn is the efficiency parameter of 
Eq. IIUI ) With this change of cooling flow accretion model 
and taking rj = 0.03, we are able to recover the results 
of Sections |3 and |1| without changing any other aspects 
of our galaxy and black hole formation model. The final 
galaxy population is, in fact, almost identical to that pre- 
sented in previous sections. This is not surprising, perhaps, 
since Eq. I29l has very similar scaling properties to Eq. 1101 In 
Fig. ll2l we illustrate the success of the model by overplotting 
its prediction for the K-band luminosity function (the solid 
line) on the observational data and on the prediction of the 
cold cloud accretion model of the last subsection. The two 
models agree very closely both with each other and with our 
standard phenomenological model (see Fig.|Sll. 



6 CONCLUSIONS 

AGN feedback is an important but relatively little explored 
element in the co-evolution of galaxies and the supermassive 
black holes at their centres. In this paper we set up machin- 
ery to study this co-evolution in unprecedented detail using 
the very large Millennium Run, a high-resolution simula- 
tion of the growth of structure in a representative region of 
the concordance ACDM cosmology. Most of our modelling 
follows earlier work, but in an important extension we in- 
troduce a 'radio' feedback mode, based on simple physical 
models and on the observed phenomenology of radio sources 
in cooling fiows. This mode suppresses gas condensation at 
the centres of massive haloes without requiring the forma- 
tion of new stars. Our modelling produces large catalogues 
of galaxies and supermassive black holes which can be used 
to address a very wide range of issues concerning the evo- 
lution and clustering of galaxies and AGN. Some cluste ring 
results were already presented in ISprineel et aI](l2005b^ . In 
the present paper, however, we limit ourselves to present- 
ing the model in some detail and to investigating the quite 
dramatic effects which 'radio mode' feedback can have on 
the luminosities and colours of massive galaxies. Our main 
results can be summarised as follows: 

(i) We study the amount of gas supplied to galaxies in each 
of the two gas infall modes discussed by IWhite fc FrenkI 
(1991): the 'static halo' mode where postshock cooling is 
slow and a quasistatic hot atmosphere forms behind the 
accretion shock; and the 'rapid cooling' mode where the 
accretion shock is radiative and no such atmosphere is 
presen t. W c distinguish these two modes using th e crite- 
rion oflWh itc & Frcnk (1991) as modified by Spring el et alJ 
('2001air and tested explicitly using SPH simulations by 
iYoshida et al.l Ill2002.'l . Our results show a sharp transition 
between the two regimes at a halo mass of 2-3 x IO'^^Mq. 
This division depends on the chemical enrichment prescrip- 
tion adopted and moves from higher to lower Vvir with time 
(corresponding to approximately constant Mvir), suggesting 
a 'down-sizing' of star formation activity as the bulk of the 
gas accreted by the haloes of larger systems is no longer 
available in the interstellar medium of the central galaxy. 

(ii) We have built a detailed model for cooling, star for- 
mation, supernova feedback, galaxy mergers and metal en- 
richm e nt based on th e earlier mod els of [ Kauffma m^^lJ 
Jl999h . ISprineel et all J2001al) and iDe Lucia etld] J2004) . 
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Applied to the Millennium Run this model reproduces many 
of the observed properties of the local galaxy population: 
the TuUy-Fisher, cold gas fraction/stellar mass and cold gas 
metallicity/stellar mass relations for Sb/c spirals (Fig. 
the field galaxy luminosity functions (Fig. |H| & IllH : the 
colour- magnitude distribution of galaxies (Fig. 0; and the 
increase in mean stellar age with galaxy mass (Fig. llOt . In 
addition the model produces a global star formation his- 
tory in reaso nable agreement with o bservation (Fig.l^. We 
also show in ISprineel et all ll2005d) that the z = cluster- 
ing properties of this population are in good agreement with 
observations. 

(iii) Ou r black hole implementa t ion e xtends the previous 
work of iKauffmann fc Haehneill l)200rl) by assuming three 
modes of AGN growth: merger-driven accretion of cold disk 
gas in a 'quasar mode', merging between black holes, and 
'radio mode' accretion when a massive black hole finds itself 
at the centre of a static hot gas halo. The 'quasar mode' is 
the dominant source for new black hole mass and is most 
active between redshifts of four and two. The 'radio mode' 
grows in overall importance until z = and is responsible for 
the feedback which shuts off the gas supply in cooling flows. 
This model reproduces the black hole mass/bulge mass re- 
lation observed in local galaxies (Fig.|lJ. The global history 
of accretion in the 'quasar mode' is qualitatively consistent 
with the evolution inferred from the optical AGN population 
(Fig.H. 

(iv) Although the overall accretion rate is low, 'radio mode' 
outflows can efficiently suppress condensation in massive 
systems (Fig.|7|l. As noted by many authors who have stud- 
ied the problem in more detail than we do, this provides an 
energetically feasible solution to the long-standing cooling 
flow 'problem'. Our analysis shows that the resulting sup- 
pression of gas condensation and star formation can produce 
luminosity functions with very similar bright end cut-OA's to 
those observed (Fig. (HJ, as well as colour-magnitude dis- 
tributions in which the most massive galaxies are red, old 
and elliptical, rather than blue, young and disk-dominated 
(Figs 1^ and ITUl. 

(v) The B — V colour distribution of galaxies is bimodal at 
all galaxy masses. Galaxies with early-type bulge-to-disk ra- 
tios are confined to the red sequence, as are the most mas- 
sive galaxies, and the most massive galaxies are almost all 
bulge-dominated, as observed in the real universe (Fig. [nj. 
This bimodality provides a natural division of model galax- 
ies into red and blue subpopulations. The colour-dependent 
luminosity functions are qualitatively similar to those found 
for early and late- type galaxies in the 2dFGRS (Fig. Illll . 
although there are significant discrepancies. After exhaust- 
ing their cold gas, massive central galaxies grow on the 
red sequence through 'burstless' merging, gaining a factor 
of two or three in mass withou t significant star formation 
iAragon-Salamanca et ai]ll998fl . Such hierarchical growth 
does not change a galaxy's colour significantly, moving it 
brightward almost parallel to the colour-magnitude relation. 

(vi) We present two physical models for black hole accretion 
from cooling fiow atmospheres. We suppose that this accre- 
tion is responsible for powering the radio outflows seen at 



the centre of almost all real cooling flows. The models differ 
in their assumptions about how gas accretes from the inner 
regions of the cooling flow, where it is thermally unstable 
and dynamically collapsing. One assumes accretion of cold 
gas clouds if these come within the sphere of influence of 
the black hole, while the other assumes Bondi-like accretion 
from the residual diffuse hot gas component. Each of the two 
models can produce z = galaxy populations similar both to 
that of our simple phenomenological model for 'radio mode' 
feedback and to the observed population (see Fig. I12II . Our 
main results are thus not sensitive to the details of the as- 
sumed accretion models. 

The presence of heating from a central AGN has long 
been suspected as the explanation for the apparent lack of 
gas condensation in cluster cooling flows. We have shown 
that including a simple treatment of this process in galaxy 
formation models not only 'solves' the cooling flow problem, 
but also dramatically affects the properties of massive galax- 
ies, inducing a cut-off similar to that observed at the bright 
end of the galaxy luminosity function, and bringing colours, 
morphologies and stellar ages into much better agreement 
with observation than is the case for models without such 
feedback. We will extend the work presented here in a com- 
panion paper, where we investigate the growth of supermas- 
sive black holes and the related AGN activity as a function 
of host galaxy properties out to high redshift. The cata- 
logues of galaxies and supermasive black holes produced by 
our modelling machinery are also being used for a very wide 
range of projects related to understanding formation, evo- 
lution and clustering processes, as well as for interpreting 
observational samples. 
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